#pragma once
#include <fftw3.h>

#include <math.h>
#include "complex_types.h"

#ifdef __CUDACC__
# include "hdev.h"
# include <cufft.h>
#endif

static inline int pad_index (int M2, int M){
	return (int) floor (M2/2.0) + (int) ceil (-M/2.0);
}


template <typename xtype> 
inline void _hdev_ ifftshift (xtype *x, int n, int s = 1)
{
	if ((n % 2) == 0)
	{
		for (int i = 0, p = (n/2); p < n; ++i, ++p)
		{
			xtype ti = x[s*i], tp = x[s*p];

			x[s*i] = tp;
			x[s*p] = ti;
		}
	}
	else
	{
		int p = (n/2);
		xtype te = x[s*(n-1)];
		xtype tp = x[s*p];

		for (int i = 0; p < n-1; ++p, ++i){
			xtype ti = x[s*i];
			x[s*i] = tp;
			tp = x[s*(p+1)];
			x[s*(p+1)] = ti;
		}

		x[n/2] = te;
	}
}

template <typename xtype> 
inline void _hdev_ fftshift (xtype *x, int n, int s = 1)
{
	if ((n % 2) == 0)
	{
		for (int i = 0, p = (n+1)/2; p < n; ++p, ++i){
			xtype ti = x[s*i], tp = x[s*p];
			x[s*i] = tp;
			x[s*p] = ti;
		}
	}
	else
	{
		int i = (n+1)/2-1;
		xtype ti = x[s*i];

		for (int p = n-1; i > 0; --i, --p){
			xtype tp = x[s*p];
			x[s*p] = ti;
			ti = x[s*(i-1)];
			x[s*(i-1)] = tp;
		}

		x[s*((n+1)/2-1)] = ti;
	}
}

template <typename xtype>
static inline void ifft2shift (xtype *x, int xs1, int xs2)
{
	for (int i2 = 0; i2 < xs2; ++i2)
		ifftshift(&x [i2*xs1], xs1, 1);

	for (int i1 = 0; i1 < xs1; ++i1)
		ifftshift(&x[i1], xs2, xs1);
}

template <typename xtype>
static inline void fft2shift (xtype *x, int xs1, int xs2)
{
	for (int i2 = 0; i2 < xs2; ++i2)
		fftshift (&x [i2*xs1], xs1, 1);

	for (int i1 = 0; i1 < xs1; ++i1)
		fftshift (&x [i1], xs2, xs1);
}

void plan_fft2c (fftwf_plan &forward, fftwf_plan &backward,
                 cplxf *X, int xs1, int xs2, int nch);

void destroy_fft2c (fftwf_plan &forward, fftwf_plan &backward);

void fft2c_exec (fftwf_plan& direction, float scale,
                 cplxf *X, int xs1, int xs2, int nch);

#ifdef __CUDACC__

void plan_cufft2c (cufftHandle &plan, int xs1, int xs2, int batch);

void destroy_cufft2c (cufftHandle &plan);

void cufft2c_exec (cufftHandle& plan, int dir,
                   cplxf *X, int xs1, int xs2, int batch);

void cuifft2shift (cplxf *X, int xs1, int xs2, int batch);
void cufft2shift (cplxf *X, int xs1, int xs2, int batch);

#endif

void ifft3c (cplxf *data, int *dsize);

